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Land surface albedo has been recognized by the Global Terrestrial Observing System (GTOS) as an essential climate 
variable crucial for accurate modeling and monitoring of the Earth's radiative budget. While global climate studies 
can leverage albedo datasets from MODIS, VIIRS, and other coarse-resolution sensors, many applications in hetero- 
geneous environments can benefit from higher-resolution albedo products derived from Landsat. We previously 
developed a “MODIS-concurrent” approach for the 30-meter albedo estimation which relied on combining post- 
2000 Landsat data with MODIS Bidirectional Reflectance Distribution Function (BRDF) information. Flere we 
present a “pre-MODIS era” approach to extend 30-m surface albedo generation in time back to the 1980s, through 
an a priori anisotropy Look-Up Table (LUT) built up from the high quality MCD43A BRDF estimates over represen- 
tative homogenous regions. Each entry in the LUT reflects a unique combination of land cover, seasonality, terrain 
information, disturbance age and type, and Landsat optical spectral bands. An initial conceptual LUT was created 
for the Pacific Northwest (PNW) of the United States and provides BRDF shapes estimated from MODIS observa- 
tions for undisturbed and disturbed surface types (including recovery trajectories of burned areas and non-fire dis- 
turbances). By accepting the assumption of a generally invariant BRDF shape for similar land surface structures as a 
priori information, spectral white-sky and black-sky albedos are derived through albedo-to-nadir reflectance ratios 
as a bridge between the Landsat and MODIS scale. A further narrow-to-broadband conversion based on radiative 
transfer simulations is adopted to produce broadband albedos at visible, near infrared, and shortwave regimes. We 
evaluate the accuracy of resultant Landsat albedo using available field measurements at forested AmeriFlux sta- 
tions in the PNW region, and examine the consistency of the surface albedo generated by this approach respective- 
ly with that from the “concurrent” approach and the coincident MODIS operational surface albedo products. Using 
the tower measurements as reference, the derived Landsat 30-m snow-free shortwave broadband albedo yields an 
absolute accuracy of 0.02 with a root mean square error less than 0.016 and a bias of no more than 0.007. A 
further cross-comparison over individual scenes shows that the retrieved white sky shortwave albedo from the 
“pre-MODIS era” LUT approach is highly consistent (R 2 = 0.988, the scene-averaged low RMSE = 0.009 and 
bias = — 0.005) with that generated by the earlier “concurrent” approach. The Landsat albedo also exhibits 
more detailed landscape texture and a wider dynamic range of albedo values than the coincident 500-m MODIS 
operational products (MCD43A3), especially in the heterogeneous regions. Collectively, the “pre-MODIS” LUT 
and “concurrent” approaches provide a practical way to retrieve long-term Landsat albedo from the historic 
Landsat archives as far back as the 1980s, as well as the current Landsat-8 mission, and thus support investigations 
into the evolution of the albedo of terrestrial biomes at fine resolution. 

© 2014 Elsevier Inc. All rights reserved. 


1. Introduction 

Surface albedo, defined as the ratio of radiant flux reflected from the 
Earth's surface to the incident flux, has been documented by the Global 
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Terrestrial Observing System (GTOS) as one of the essential climate var- 
iables governing Earth's surface energy budget (Pinty et al., 2008; 
Schaaf, Cihlar, Belward, Dutton, & Verstraete, 2009; Schaaf et al, 
2008). The radiative forcing intercepted by the land surface is perhaps 
the most important initial energy source for biophysical processes, 
through a further conversion into latent, sensible, and stored heat 
terms and input to the soil-vegetation biophysical system (Betts, 
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2000; Lyone, Jin, & Randerson, 2008; Ollinger et al., 2008; Peckham, Ahl, 
Serbin, & Gower, 2008; Randerson et al, 2006; Sellers, Los, et al., 1996; 
Sellers, Randall, et al., 1996; Zhang et al., 2009). Studies have shown 
that land cover change (and ecosystem disturbance) may have a signifi- 
cant influence on regional albedo, and hence long-term climate forcing 
(Bala et al., 2007; Betts, 2000; Claussen, Brovkin, & Ganopolski, 2001; 
Randerson et al., 2006). Terrestrial albedo varies enormously in space 
and time as a result of both natural events (e.g. weather disaster, insect, 
disease, wild fire, season-shifts, and vegetation phenological phase) and 
human activities (e.g. forest-thinning & clearing, crops-sowing & harvest- 
ing, urbanization, and other land use management methods) (Jin & Roy, 
2005; Ju, Roy, Shuai, & Schaaf, 2010; OfHalloran, et al., 2011; Shuai and 
Schaaf, 2010; Shuai, Schaaf, et al., 2013; Shuai, Xie, Wang, & Wu, 2013; 
Xu et al., 2013). As strategies emerge for managing ecosystem carbon in 
order to mitigate global warming, several studies have pointed out the 
potential risk of ignoring the physical consequences of land cover change, 
including changes to land surface albedo (Betts, 2000; Lyone et al., 2008; 
Peckham et al., 2008; Randerson et al., 2006). 

Albedo datasets have been derived from existing coarse-resolution satellite 
sensors to parameterize global land surface and climate models. Compared 
with previous single-angle models, modem albedo algorithms rely on multiple 
directional reflectance measurements to first estimate a Bi-directional Reflec- 
tance Distribution Function (BRDF) model of the target, then integrate over in- 
cident and view hemispheres to calculate albedo. Studies have concluded that 
relative errors can reach up to 45% without the consideration of direction/ 
angle effects in the albedo estimation (Kimes & Sellers, 1985; Kimes, Sellers, 
& Newcomb, 1987). Because most satellite sensors cannot collect multiple ob- 
servations of a target in a single pass, the sequential accumulation of data over 
multiple days (for sun-synchronous orbit) or multiple hours (geostationary 
orbit), may be adopted as a relevant solution to achieve multi-angle measure- 
ments sampling the full sun-target-sensor geometry. Global surface albedo 
has been mapped from the Advanced Very High Resolution Radiometer 
(AVHRR) (Csiszar & Gutman, 1999; Key, Wang, Stroeve, & Fowler, 2001), 
Earth Radiation Budget Experiment (ERBE) radiometer data (Li & Garand, 
1994), and the Along Track Scanning Radiometer (ATSR). With the advent 
of routine albedo products retrieved from Polarization and Directionality of 
the Earth's Reflectances (POLDER-I and II) (Bicheron & Leroy, 2000; 
Hautecoeur & Leroy, 1998; Leroy et al., 1997; Maignan, Breon, & Lacaze, 

2004) , Multi-angle Imaging SpectroRadiomenter (MISR) (Martonchik, Pinty, 
& Verstraete, 2002; Martonchik et al., 1998), Clouds and the Earth's 
Radiant Energy System (CERES) (Rutan et al., 2009), Meteosat Visi- 
ble and Infrared Imager (MVIRI)/Meteosat and Meteosat Second 
Generation (MSG) (Carrer, Roujean, & Meurey, 2010; Geiger, 
Carrer, Franchisteguy, Roujean, & Meurey, 2008; Pinty et al., 2000), 
SPOT4/VEGETATION (Franchisteguy, Geiger, Roujean, & Samain, 

2005) , and the recently launched Visible Infrared Imager Radiometer 
Suite (VIIRS) (Justice et al., 2013; Liang, Yu, & Defelice, 2005), albedo 
maps with spatial resolutions of 500-m to tens of kilometer and tempo- 
ral frequencies of daily to monthly are now available to serve for climate 
model refining and inter-annual exploration (Schaaf et al., 2008). 


While global climate studies can utilize the coarse-resolution surface 
albedo datasets described above, there remains a need for consistent, 
fine-resolution albedo products for specific applications. Several publi- 
cations have highlighted the importance of land cover change, including 
deforestation, afforestation, agricultural expansion, urbanization, and 
other human-induced land surface alteration, to the terrestrial carbon 
cycle and climate changes (Goward et al., 2008; Masek & Collatz, 
2006; Pan et al., 201 1 ; Randerson et al., 2006). However, spatial resolu- 
tions coarser than 250-m may be insufficient to capture patch-scale 
vegetation changes associated with human land use and forest distur- 
bance (Townshend and Justice 1988; Masek etal., 2013). Fine resolution 
imagery (~30 m or better) can more accurately quantify the areas and 
rates of these anthropogenic land changes. In addition, for climate 
change investigations, long time series of albedo products are required. 
Although operational albedo datasets covering the last 30 years have 
been assembled from different sensors covering different periods, 
the merging of multiple records raises issues of data consistency and 
quality. Because of the differences among sensors (wavelength of spectral 
bands, orbit geometry, spatial resolution, and geographic region), the de- 
rived albedo products may differ depending on the specific product, the 
data source, and the production strategies (Schaaf et al., 2009). Therefore, 
datasets derived from a single continuous acquisition program offers a 
greater potential for consistency in data quality. Despite differences in 
sensor design over time, the Landsat program has acquired a 42 -year re- 
cord of Earth Observations that captured global land conditions and dy- 
namics through six successful missions since 1972. With the launch of 
Landsat-8 in February 2013 (Loveland & Dwyer, 2012), this record has 
the potential of reaching 50 years. The opening of the Landsat ar- 
chive for free distribution in late 2008 has invigorated the push for 
creating long-term biophysical and land cover products from new 
and archived Landsat data (Woodcock et al., 2008; Wulder, Masek, 
Cohen, Loveland, & Woodcock, 2012). It includes this effort to devel- 
op the long-term, consistent surface albedo products from the 
Landsat program. 

In a previous study, we developed a “concurrent” approach for gen- 
erating 30-m resolution albedo products for the post-2000 (MODIS) era 
by combining Landsat surface reflectance with MODIS surface anisotro- 
py information (Shuai, Masek, Gao, & Schaaf, 2011). In this study, we 
propose and validate a new approach to generate Landsat albedo prod- 
ucts for the pre-MODIS era, by using albedo-to-nadir reflectance ratios 
(Shuai et al., 2011) and an a priori anisotropy Look-Up Table (LUT) 
that has been built up from the high quality MCD43A BRDF retrievals 
over representative homogeneous regions. This approach yields 
both spectral and broadband albedos, and a quality assessment (QA) 
map based on the quality of MODIS anisotropy and Landsat surface re- 
flectance. In this paper, we first address the theoretical basis of the 
“pre-MODIS-era” LUT approach, creation of the BRDF-LUT, and then 
demonstrate its application over more than 100 Landsat scenes in the 
Pacific Northwest of the United States where simultaneous ground 
measurements are available for validation. 


2. Albedo definition 

The spectral Directional-Hemispherical Reflectance (DHR) of a plane surface is defined as the ratio of radiant energy scattered upward from the 
surface in all directions to the down-welling incident irradiance on the surface within the target spectrum regime ( Ai, A 2 ). It equals the integral of the 
BRDF over the view hemisphere for an incident beam at a given wavelength, as shown in formula (1). Under the extreme condition that no diffuse 
radiation but only the direct beam arrives from the solar incidence angle (0, <p) defined by zenith angle 0, and azimuth angle <p (L(0, <p)), the albedo is 
referred to as “Black-Sky Albedo” (BSA) R(0j, <p f ; A) in the MODIS product series (Lucht, Schaaf, & Strahler, 2000; Strahler et al., 1999). Under the as- 
sumption that alHrradiance is isotopic (purely diffuse skylight), a further integral over illumination hemisphere provides the Bi-Hemispherical Re- 
flectance (BHR) R( A), or “White-Sky Albedo” (WSA) formulae (2) and (3) (Lucht et al., 2000; Strahler et al., 1999). The spectral BHR under actual 
atmospheric conditions (known as the “blue-sky albedo”, or “actual albedo”) can be approximated through a linear combination of BSA and WSA, 
weighted by the fraction of actual direct to diffuse skylight (Lewis & Barnsley, 1994; Lucht et al., 2000; Romanet al., 2010). Because the upwelling 
radiance depends on not only the BRDF properties of the observed surface, but also atmospheric conditions, R( A) may change with the variation 
of the instantaneous cloud cover and aerosol loading, as well as over the course of the day as the solar geometry changes even for constant atmo- 
spheric and surface conditions (Lucht et al., 2000). In addition, multiple scattering between surface and atmosphere affects the angular distribution 
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of sky radiance. Therefore, bi-hemispheric reflectance (i.e. albedo) is not a true surface property, but rather a function of solar beam direction, 
atmospheric state, and surface anisotropic features. 


2TITT/2 

A) = J J / r (0i, <p,-; 6 V , <p v ; A) cosd v sin6 v dd v d(p v 

0 0 
2tt rr/2 

= -J J R(di,(Pj;O v ,(p v ; A) cos6 v s'md v dd v d(p v 
0 0 


( 1 ) 


where R(0 f , <p f ; A) = Spectral black-sky albedo (Directional-Hemispherical Reflectance, DHR) as a function of the solar incidence angle (0;, (pi) 
(Strahler et al., 1999), and f r (d it 0 V , c p v ; A) = Bidirectional Reflectance Distribution Function (BRDF) describing the behavior of surface scattering 

as a function of a parallel incident beam from one direction (0/, c pi) in the illuminating hemisphere into the reflected direction (0 V , <j> v ) in the viewing 
hemisphere, at a particular wavelength A. Further elaboration is presented in Nicodemus, Richmond, Ginsberg, and Limperis (1977) and Schaepman- 
Strub, Schaepman, Painter, Dangel, and Martonchik (2006). The terms “BRDF” and “anisotropy” in this paper refer to this underlying property. 
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Strictly speaking, for natural targets, BRDF or any other nominal directional-related metric is not a measurable quantity, as it requires perfectly 
collimated beams of illumination and observation, while actual sunlight is partly diffuse and the measurements involve conical geometries. Thus, in- 
dividual satellite measurement provides only an approximation of the directional reflectance. 

For most of the applications involving energy balance, the reflectance quantity of interest is not the spectral reflectance but rather reflectance in- 
tegrated over a broad spectral interval (Ai, A 2 ), see formula (4), to capture the overall radiative forcing. The spectral integrals for the hemispherical 
reflectance are functions of the down-welling solar spectrum as defined in the above formulae. The visible regime (0.3-0.7 pm) known as 
photosynthetically-active radiation (PAR) is of special interest to carbon cycle modelers for the estimation of carbon fixation via photosynthesis 
(Dorman & Sellers, 1989). In contrast, the total shortwave regime (0.3-3.0 pm), as well as visible and near-infrared bands, are typically required 
by surface energy balance studies. Note that the generic term “albedo”, without any specification of the sun-view geometry and integral wavelength, 
often implies the bi-hemispheric broadband albedo of the whole solar irradiance domain. 


R(A^A 2 ) = 
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( 4 ) 


3. Algorithm 

The initial impetus to develop a “pre-MODIS era” approach flowed 
from the desire to understand albedo consequences of specific types 
of forest disturbance and recovery at a fine resolution. To extend 30-m 
surface albedo generation in time back to the 1980s, we first build 
an a priori anisotropy Look-Up Table (LUT) from the high quality 
MCD43A BRDF estimates over representative homogenous regions, 
then calculate the albedo-to-nadir reflectance ratios for each entry 
and apply these ratios to the 30-m Landsat nadir reflectance. Finally, 
we use narrow- to broad-band conversion factors to derive broadband 
Landsat albedos. Fig. 1 outlines the overall workflow of this approach 
into three main functional components: surface reflectance calculation 
and assessment (Fig. 1A), BRDF-LUT creation (Fig. IB), and Landsat sur- 
face albedo generation (Fig. 1C). The aim of “surface reflectance calcula- 
tion and assessment” is to retrieve terrain and atmosphere corrected 
surface reflectance (that is defined in Masek et al., 2006) from the 


Landsat Thematic Mapper and Enhanced Thematic Mapper Plus 
(TM/ETM + ) LIT, and remove pixels contaminated by cloud and 
snow from further analysis. The aim of “BRDF-LUTs creation” is to 
build up the a priori anisotropy information for each defined land 
surface category from the operational MODIS 500-m high quality an- 
isotropy products (i.e. MCD43A) over representative homogeneous 
land surface structure regions. The aim of “Landsat albedo generation” 
is to obtain the narrow-band spectral albedo by combining Landsat di- 
rectional surface reflectance with the specific a priori anisotropy infor- 
mation stored in the BRDF-LUTs, and then convert narrow to broad 
band albedos for the visible (0.3-0.7 pm), NIR (0.7-3.0 pm), and short- 
wave (0.3-3.0 pm) regimes. 

3.1. Surface reflectance assessment 

Landsat surface directional reflectance at each spectral band of The- 
matic Mapper and Enhanced Thematic Mapper Plus (TM/ETM + ) has 
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Fig. 1 . Flow chart of the “pre-MODIS era” LUT approach composed of three functional components (A. surface reflectance assessment; B. BRDF-LUT creation; and C. Landsat albedo 
generation). 


been produced from orthorectified Landsat level 1 T raw images 
downloaded from USGS EROS , using the Landsat Ecosystem Disturbance 
Adaptive Processing System (LEDAPS) (Masek et al M 2006). The level 1 
raw radiometry data ingested by LEDAPS were calibrated to at-sensor 
radiance, converted to top-of-atmosphere reflectance, and then atmo- 
spherically corrected to surface reflectance using the second simulation 
of satellite signal in the solar spectrum (6S) model (Masek et al., 2006; 
Vermote, Saleous, & Justice, 2002; Vermote et al., 1997). LEDAPS dem- 
onstrated good performance through comparisons with ground-based 
AERONET optical thickness measurements (Masek et al., 2006), concur- 
rent MODIS Terra reflectance (Feng et al., 2012; Masek et al., 2006), and 
other approaches for Landsat surface reflectance generation (Ju, Roy, 
Vermote, Masek, & Kovalskyy, 2012). To mitigate the cloud effect on 
the surface radiometric accuracy, pixels contaminated by cloud, cloud 
shadow, and adjacent clouds were screened from this study using the 
LEDAPS-derived cloud mask. An additional screening for snow was per- 
formed based on the operational MODIS snow mapping algorithm (Hall, 
Riggs, Salomonson, DiGirolamo, & Bayr, 2002), through the Normalized 
Difference Snow Index (NDSI) calculated from reflectance at Landsat 
green (0.53-0.61 pm) and shortwave infra-red (1.55-1.75 pm) bands. 
Further thresholds for green band reflectance (>0.10) and NDVI were 
applied to reduce the erroneous classification of very dark targets 
(such as black spruce forests), as well as the thermal mask to eliminate 
the spurious snow cover possibly induced by residual cloud cover, aero- 
sol effect and snow/sand confusion on coastlines (Hall, Riggs, & 
Salomonson, 1995; Hall et al., 2002). 

3.2. BRDF Look-Up Tables 

The most direct way to obtain anisotropy information of any land 
surface target at the pixel scale is to collect a representative sample of 
reflectance observations at multiple directions, over a short interval of 
time. However, because of the narrow field of view of Landsat ( ± 7.5 de- 
grees) and the limited number of acquisitions offered by the 16-day re- 
peat cycle, it is not feasible to obtain target anisotropy information 


directly from multiple Landsat directional reflectance observations. In- 
stead, we need to obtain target BRDF estimates from other sources, 
such as MODIS, or MISR. 

For this study, the Collection V005 MODIS 8-day anisotropy dataset 
(MCD43A) was used to create the BRDF-LUT because of its wide range 
of sun and view angles, the broad spectral coverage of MODIS for simul- 
taneous atmosphere correction, frequent acquisition for the potential 
daily adjustment of BRDF retrieval, the 500-m moderate resolution, 
and especially the continuity of global products since 2000. The opera- 
tional MODIS albedo and reflectance anisotropy products make use of 
the kernel-driven, linear algorithm that relies on the weighted sum of 
an isotropic and two additional kernels (respectively called Ross-thick 
and Li-sparse-reciprocal models, RTLSR) of viewing and illumination ge- 
ometry to estimate the BRDF model (Li & Strahler, 1992; Lucht et al., 
2000; Ross, 1981 ; Roujean, Leroy, Podaire, & Deschamps, 1992). The re- 
trieved kernel weights (also called BRDF model parameters) are those 
that best fit an adequate angular sample of the high quality cloud- 
cleared, atmospherically corrected surface reflectances available for 
each pixel over a 16-day period (Lucht et al., 2000; Schaaf et al., 
2002, Schaaf, Liu, Gao, & Strahler, 2011; Shuai, Schaaf, Strahler, Liu, 
& Jiao, 2008; Shuai & Schaaf, 2010). This model combination has 
been shown to be well-suited to describe the surface anisotropy of 
the variety of land surfaces distributed worldwide (Privette, Eck, & 
Deering, 1997). The absolute accuracy of MCD34A albedo at local 
solar noon (LSN) derived from the estimated BRDF model has been 
established by comparison with ground measurements from avail- 
able international Baseline Surface Radiation Network (BSRN) and 
Fluxnet sites (Cescatti et al., 2012; Roman et al., 2009; Wang et al., 
2014). This algorithm assumes that the land surface does not experi- 
ence significant structural changes during the 16-day observation 
period, which is reasonable except in circumstances of abrupt distur- 
bance or conversion. 

The creation of a BRDF LUT is based on the identification of land sur- 
face intrinsic anisotropic features which make one object distinguishable 
from others. Numerous studies have demonstrated unique anisotropic 
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Table 1 

Structure of the BRDF LUT. 


Type 

Land cover class # 

Disturbance age 

Disturbance severity 

Month 

DEM 

QA 

Bands (1-5,7) 

RTLS-R parameters 

Range of value 

NLCD classification scheme 

0-30 

Low, medium, or high 

1-12 

Flat or mountainous 

0-5 

Landsat 1-5, 7 

Isotropic, volumetric, and 
geometric kernel weights 

Un-disturbed 


NA 

NA 






Fire-disturbed 

\S 

\S 



NA a 


\S 


Nonfire-disturbed 



NA a 


NA a 


\S 



a Limited by the lack of current ancillary data community. 


features among distinct landscape attributes (Bacour & Breon, 2005; 
Bicheron & Leroy, 2000; Lovell & Graetz, 2002; Maignan et al M 2004; 
Shuai & Schaaf, 2010; Strugnell & Lucht, 2001), biome components 
(Chen & Leblanc, 1997), vegetation life-cycle and seasonal stages 
(Kimes, 1983; Shuai et al., 2011), and classes of disturbance induced by 
natural or human activity, as well as significant effects from terrain 
(Schaaf, Li, & Strahler, 1994). Thus, the attributes defined in Table 1 are 
adopted to build the overall conceptual structure of BRDF-LUT. Each 
entry in the LUT reflects a unique combination of land surface type, ter- 
rain, time of year, limited disturbance age and type, and Landsat spectral 
bands. 

Several ancillary datasets provided the basis for this stratification. 
First, the 30-m 2006 NLCD (National Land Cover Database, Vogelmann, 
Sohl, & Howard, 1998) classification maps with high overall and user's 
accuracy (Wickham et al., 2013) were used to determine local landscape 
attributes, and to identify representative homogenous land surface re- 
gions when aggregated to the MODIS 500-m resolution. Then, two 
datasets giving the timing and location of ecosystem disturbance were 
used to quantify the BRDF evolution of disturbed landscapes. The annual 
30-m Monitoring Trends in Burn Severity (MTBS) (Eidenshink et al., 
2007) dataset has mapped the low/medium/high burn severity of 
fires (greater than 1000 acres in the west and 500 acres in the east) 
that have occurred since 1984 across all lands of the United States. The 
30-m NAFD (North American Forest Dynamics, Masek et al., 2008; 
Masek et al., 2013; Huang et al., 2010) dataset identified other forest 
non-fire disturbance events (such as harvest, storm damage, or disease) 
over the same time period. While the NAFD dataset targets rapid distur- 
bance events that remove substantial canopy cover, more subtle or grad- 
ual declines in live biomass (e.g. selective tree removal, gradual insect 
outbreaks) may not be captured. Both MTBS and NAFD datasets are gen- 
erated from Landsat spectral signatures before and after the disturbance 
events. While the MTBS dataset uses independent confirmation of fire 



VZN °( degrees) 


Fig. 2. Example of the difference in MODIS BRDF shape estimated for non-disturbed ever- 
green forest (in principle plane at 30° solar zenith angle) obtained from a mountainous re- 
gion (slope > 15°, dot-line) and a relatively flat region (slope < 15°, solid line) at NIR 
(upper), SWIR (middle), and Red (lower) bands from the Pacific Northwest region of 
the United States. 


timing (from the National Interagency Fire Center), the NAFD dataset 
may have one or two year bias on the timing of disturbance if a cloud- 
free image or composite was not produced for a given year. Since both 
datasets cover the period since 1984, any fire and non-fire disturbance 
events encountered before the Landsat TM/ETM + era (1980s) are not 
be defined in the maps. In addition, the Shuttle Radar Topography Mis- 
sion (SRTM) DEM data (Farr et al., 2007) was utilized to differentiate the 
quality of BRDF shapes affected by mountainous terrain. Due to the 
study in Schaaf et al. (1994) in terrain effects on the anisotropy feature, 
the MODIS estimated BRDFs were qualitatively graded into two strata 
(flat to moderate with slope <15°; and steeper mountainous with 
slope >15°). An example is shown in Fig. 2 for individual undisturbed 
evergreen needle leaf forest patches respectively in a flat region and a 
high-slope mountainous region. 

Thus the a priori BRDF LUTs for the Pacific Northwest (PNW) region 
of United States were created from MCD43A products and ancillary 
datasets (see module B in Fig. 1 ) in terms of the above conceptual struc- 
ture of the BRDF LUT. The PNW was selected for this initial prototype 
due to its range of ecosystems, prevalence of both fire- and non-fire 
forest disturbances, and range of topography. In order to minimize the 
effect of biome mixtures, MODIS 500-m pixels were labeled as repre- 
sentative “pure” pixels if they were composed of at least 85% of a single 
land surface type when aggregated from the 30-meter NLCD land cover 
map. In addition to being stratified by NLCD land cover, the LUT of BRDF 
was also stratified by disturbance type (“undisturbed”, “fire disturbed”, 
“non-fire disturbance”), disturbance severity (from the MTBS fire dis- 
turbance product), topographic slope (greater or less than 15°), time 
since disturbance (0-26 years corresponding to the NAFD and MTBS 
coverage of 1984-2010), and month of the year (Table 1 ). For each com- 
bination of these attributes, the BRDF shapes for Landsat (and MODIS) 
reflective bands were extracted from the operational V005 8-day 
MCD43A1 (BRDF parameters) and MCD43A2 (QA flags) 11 -year prod- 
uct (Schaaf et al., 2002; Shuai et al., 2008). The time dimension 
(month for the undisturbed LUT, and age of disturbance), was used to 
depict the seasonality, growth phase, and growth evolution since distur- 
bance, in the BRDF shapes over for a given land surface scenario. If no 
high quality BRDF was available for a given month (for a seasonal char- 
acterization) or year (for characterizing post-disturbance evolution), a 
backup BRDF shape was established through linear interpolation of 
the BRDF model parameters from available time periods. To document 
the quality of BRDF shapes in the LUT, each was assigned a quality flag 
denoted as “high quality” for the original MCD43A estimation and 
“low quality” for those interpolated ones. 

As an example, Fig. 3 shows the BRDF shapes in the principle plane 
with solar incident at 30° zenith angle, averaged over disturbed ever- 
green forest regions in the PNW. The snow-free time series of BRDF 
shapes from September illustrate the evolution of evergreen forest sig- 
nature over two decades in the green, NIR, and SWIR bands. It is seen 
that BRDFs of both fire and non-fire disturbance types have systematic 
temporal variations in shapes and magnitudes. The evolution of this 
generalized BRDF-shape may be associated with regrowth and recovery 
of canopy greenness and structure for the disturbed forest land, indicat- 
ed by the gradual sharpening or flattening of the hot-spot. There are 
strong temporal signatures of green vegetation in both examples of 
the disturbance types, firstly displayed as a clear enhanced hot-spot in 
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Fig. 3. Two decades of BRDF evolution following non-fire disturbance (harvest, thinning dominated, top panels) and high-severity fire disturbance (bottom panels) in the Pacific Northwest 
of the United States, for green (left), near-infrared (middle), and shortwave-infrared (right) bands. The original MCD43A BRDF shapes were retrieved from September in the principle 
plane with solar incident at 30° zenith angle. The BRDF shapes show a strong hot-spot in the backward (showing as positive view zenith angle— VZA) direction, and systematic changes 
in magnitude and shape following disturbance events. 


the NIR band with the increasing of greenness, and a suppression of 
the hot-spot at the SWIR with the augmentation of canopy water 
content accompanying the forest regrowth. In contrast to the mono- 
tonic decrease in brightness following non-fire disturbance, the 
high-burn-severity fire disturbed forest presents a complicated trajec- 
tory of anisotropy development in the green band, with a multiple- 
year (-7-12 years since disturbance) reduction during the increase of 
the hot-spot. It may be explained by the different recovery trajectories 
of the post-fire residual structures. These post-fire residual transition 
rates will vary among fires, with a high rate in the first two years from 
tree to snag (i.e. tree mortality), and a late peak after several years 
later for the tree-to-downed wood and snag-to-downed wood change 
depending on the species and tree size of the burned forest region. 
Once the green signature from the re-grown forest and understory veg- 
etation (such as grass or shrub) becomes dominant, a continuous grad- 
ual increasing can be captured generally 10 years after severe fires, as 
shown in Fig. 3. Some small fluctuations found in the gradual evolution 
of each BRDF shape could be due to uncertainties in the mapped timing 
of disturbance, poorer quality BRDF estimation, variations in atmo- 
spheric conditions, and residual cloud and snow effects. 


Table 2 

Segments of the pixel-based 16-bit QA word for each Landsat albedo map to indicate the 
performance of albedo retrieval. 


Bit 

Meaning 


bl5 

Fill value (1 = fill- value) 


bl4 

Cloud flag (1 = cloud contamination) 


bl3 

Snow flag (1 = snow contamination) 


bl2 

Disturbance flag (0 = undisturbed; 1 = disturbed) 


bl 1-10 

Disturbance type (00 = fire; 01 = non-fire; 10 and 11 = 

reserved) 

b9-8 

Fire disturbance severity (00 = reserved; 01 = low; 10 = 
11 = high) 

medium; 

b7 

BRDF QA (0 = original; 1 = backup/interpolation) 


b6-0 

Disturbance age for disturbed pixel or land cover class for 
un-disturbed pixel 



3.3. Surface albedo determination 


Once the BRDF shape is determined, surface albedo at 30-meter 
resolution can be calculated from the albedo-to-nadir reflectance ratio 
(A/N) and Landsat surface reflectance as detailed for the “concurrent” 
approach in Shuai et al. (2011 ). This method assumes that a given sur- 
face type has the same BRDF shape at MODIS or Landsat resolution, 
and can be scaled to albedo using the 30-meter directional reflectance 
from Landsat as shown in (5), with R /nd and R mod denoting the corre- 
sponding spectral reflectance from Landsat and MODIS, respectively. 
Then, the Landsat black-sky albedo with a solar zenith angle at the 
Landsat overpass time and white-sky albedo were computed respec- 
tively for the six non-thermal Landsat bands. The broadband albedos 
for visible (0.3-0.7 pm) a vis , , near infrared (0.7-3.0 pm) a niri and short- 
wave (0.3-3.0 pm) ois hort bands were produced by a further conversion 
from narrow spectral band albedo values (a/) using new conversion co- 
efficients for Landsat 5 TM (6-8) and Landsat 7 ETM +(9-11). These co- 
efficients were derived from radiative transfer simulations using 245 
surface spectra representing different surface types (He, Liang, Wang, 
Shuai, & Yu, 2013; Liang, 2000). Finally, a quality assessment (QA) 
layer constructed into a 16-bit word was stored for each pixel 
(see Table 2) to track the quality of input data, and estimate error prop- 
agation through the fusion of multiple data sources. 


f ^mod (®mod ~ ®*Pmod ~ *Pi 5 j *Pv 5 'M ^R-lnd i^lnd ~ j ^Plnd ~ ^Pi 5 ? *Pv ■> 'M 

I fr_ m J0i, cfc; e v , <p v ; \) = f r _J6 h cfc; 0 V , c/> v ; A) 



Rlnd(0i,<Pi-A,<Pv) >< 


Rmod(^) 

Rmod^vViAiVv) 


Rlnd^vViM 


Rlnd(Pi><Pi m ,0v,<Pv) x 


Rmod(0v<Pf0viVv) 


( 5 ) 


a shon = 0.3206oq + 0.1572a 3 + 0.3666ct 4 + 0.1162a 5 
+ 0. 045 7a 7 — 0.0063 


( 6 ) 
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Fig. 4. Examples of the “pre-MODIS era” approach generated from scene (path/row: 45/29) on day 2007-08-29. (A) The spectral black-sky albedo composite of Landsat-5 bands 5, 4, and 3, 
(B) the broadband black-sky albedo composite of visible, near infrared, and shortwave bands, (C) the black-sky albedo for the shortwave band, and (D) the quality assessment maps. 


a vis = 0.6000^ + 0. 2204a 2 + 0.1828a 3 -0.0033 (7) 

a nir = 0.6646a4 + 0.2859a s + 0.0566a 7 — 0.0037 (8) 

otshort = 0.3141a-! + 0.1607ct 3 + 0.3694a 4 + 0.1160a 5 

+ 0.0456a 7 -0.0057 (9) 

a vis = 0.5610^ + 0. 2404a 2 + 0.201 2a 3 -0.0026 (10) 

a nir = 0.66680:4 + 0.2861a5 + 0.0572a 7 -0.0042 (11) 


3.4. Central Oregon example for the derived Albedo and QA maps 

Fig. 4 shows maps of the 30-m Landsat albedo products generated 
from the “pre-MODIS era” LUT approach for a scene in central Oregon 
(path/row: 45/29) on August 29, 2007. Spectral black-sky albedo 
estimates are provided as the composite of shortwave infrared, near in- 
frared, and red bands (wavelength centered 1.65 pm, 0.83 pm, and 
0.66 pm) (Fig. 4A). Broadband black-sky albedos are available for 
the visible (0.3-0.7 pm), near-infrared (0.7-3.0 pm), and shortwave 
(0.3-3.0 pm) bands (Fig. 4B and C). At the date corresponding to the se- 
lected sample case, a large part of the Central and Eastern region was 
dominated by sparse shrubs or barren land. Compared with the forest 
region in the central-west part, these areas appear as high values in 
the SWIR and Red bands, lower values in the NIR band, with scattered 


agricultural fields in the lower-central region and forest stands in the 
middle-eastern region in Fig. 4A. In the three-broadband composite 
image (Fig. 4B), however, the albedo in the visible regime has higher 
values than the other two bands and shows up as brown-red in the cor- 
responding areas. For the retrieval of each pixel, one corresponding QA 
map (Fig. 4D) provides the possible cloud and snow contamination, and 
details of undisturbed or disturbed information. 

4. Accuracy assessment of the Landsat albedo products 

Three approaches have been used to evaluate the accuracy of albedo 
products generated by the “pre-MODIS era” LUT approach presented in 
this paper. One is the direct validation of shortwave albedo with actual 
ground measurements. The other two methods are cross-comparisons 
of surface albedo maps generated by (1) the “concurrent” approach of 
Shuai et al. (2011) that uses coincident MODIS products to retrieve 
Landsat-scale albedo, and (2) the coincident operational MODIS albedo 
products themselves. Comparison with ground measurements is an inde- 
pendent and optimal approach for product validation, but suffers from 
the limited availability of ground albedo-meter measurements. Cross- 
comparison with other products can be performed on a large volume of 
MODIS images, but does not provide a robust estimate of absolute accura- 
cy. Utilization of these multiple validation means may increase the ability 
to evaluate the algorithm performance thoroughly and objectively. 

4.1. Validation with ground measurements 

Independent ground or tower albedo measurements are generally 
considered to be more accurate than satellite retrievals, and are often 
taken as a reference for the validation of satellite products. However, 
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Table 3 

Forested ground stations in the Pacific Northwest region.Acquired from the network-wide AmeriFlux database. 


Site name 

Vegetation type 

Location 3 

Tower height (m) 

Canopy height (m) 

Footprint of 
observation (m) b 

Data period 

Land sat 
retrieval # 

US-Me2 

ENF C 

44°27 / 8.28"N, 121°33'25.92"W 

32.0 

~22 

228.6 

2005-2007; 2009-2011 

32 

US-Me3 

ENF C 

44°18'55.68"N, 121°36'28.29"W 

18.0 d 

~3.11 

342.9 

2004-2009 

4 

US-Me6 

ENF C 

44°19'23.43"N, 121°36'15.69"W 

18.6 e 

~7.0 

265.2 

2010-2011 

7 

US-NR1 

MF f 

40°1 / 58.31"N, 105°32 , 49.09"W 

26 

11.5 

331.5 

2006-201 l g 

24 

US-GLE 

Subalpine, alpine 

41°21 '59.51 "N, 106°14 / 23.82"W 

23/30 h 

12.1 

249.2/409.2 

2004-2011 

24 

US-Blk 

Conifer 

44°09 / 01"N, 103°38 / 24"W 

24 

13-15 

251.5-205.7 

2004-2009 

24 


a Location of each site is confirmed by their Pis via private communication. 
b Diameter of ground measurements footprint in the horizontal plane at canopy height. 
c Evergreen needle leaf forest. 

d Tower height is 18 m, instrument CNR-1 is mounted at 14 m. 
e Tower height is 18.6 m, instrument mounted at 17.7 m. 
f Subalpine mixed coniferous forest. 

g Currently only post 2005 ground data to be used in terms of data processor's suggestion via personal contact. 
h Tower height is 30 m during 1999-2006, and adjusted to 23 m since 2006. 


the validation of satellite-derived products is difficult because the foot- 
print of satellite observations differs significantly from that of in-situ in- 
struments. Only measurements spatially representing the surrounding 
landscape at both in-situ and satellite scales can provide a comparable 
basis for validation (Roman et al., 2009). 


4.1.1. Surface albedo ground measurements 

Tower-based surface albedo measurements were acquired from six 
available forested sites of AmeriFlux network in the Pacific Northwest 
region of the United States (Table 3; Ruehr, Martin, & Law, 2012; 
Thomas et al., 2009; Vickers, Thomas, Pettijohn, Martin, & Law, 2012; 
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Fig. 5. Distribution of Ameriflux validation sites in the Pacific Northwest region of the United States. For each site, a ground photo (Upper-left), photo of tower surroundings (lower-left), 
and high- resolution satellite image (right) are shown. Note: image of tower surroundings for the currently deactivated US-Blk site is not available. 
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Wilson & Meyers, 2007). Forested sites were of particular interest since 
one of the overall objectives was to understand how forest disturbance 
and recovery influenced the albedo trajectories. The field sites sample 
forest ecosystems with different species composition, age and distur- 
bance regimes (see the distribution and landscapes in Fig. 5), including 
sub-alpine forest older than 400 years with dispersed younger trees at 
the US-GLE site, subalpine mixed coniferous forest naturally regrown 
from extensive logging during 1900-1910 at US-NR1, conifer forest 
with scattered herbs and shrubs recovering from logging activity in 
the early 1900s at US-Me2, very young ponderosa pine stand disturbed 
by fire and harvest in the 1 980s at US-Me3, and a reforested 20-year old 
ponderosa pine site following fire and salvage cutting at US-Me6. 
Both upward and downward broadband shortwave solar radiation 
(0.3-2.8 pm) were measured via Kipp and Zonen (CNR1, CM-3, or 
CM-6b), or Eppley-PSP tower albedo-meters with 170° effective field 
of view. Data series collected from the individual sites were processed 
into the 30-minute standard values, and obtained from the Ameriflux 
web site: http://ameriflux.ornl.gov. For this study, daily tower albedo 
values were retrieved corresponding to the Landsat imaging time of 
10:30 AM, as well as local solar noon (LSN), the time corresponding to 
the MODIS MCD43A product suite. In addition, the surface albedo data 
series (level 2) for site US-NR1 were reviewed for inconsistency during 
period pre- and post-2005 as a new CNR1 sensor was installed in the 
fall of 2005. The post-2005 data which were measured with the new 
well-calibrated sensors were recommended for use by the data 
provider-(S. Burns, personal communication). An effort to establish 
sensor-to-sensor cross-calibration is underway, and may provide 
corrected pre-2005 data soon for further validation activities at US_NR1 . 

4.1.2. Aggregation from Landsat scale to tower measurement footprint 
The disparate spatial scale between satellite and in-situ measure- 
ments is one of the barriers to validating satellite-derived products. 
Several studies have concluded that direct “point-to-pixel” comparison, 
without considering spatial scales, is not sufficient for albedo product 
validation, unless the validation focuses on a large and homogenous re- 
gions (Liang et al., 2002; Roman et al., 2009). The tower based instru- 
ment pyranometer is influenced by the “cosine-law” of the response 
direction and has a 170° effective field of view. An area of 2 h x tan 
(^) diameter in the horizontal plane at forest canopy height is then de- 
fined by the downward-looking sensor mounted on a tower ( h meters 
above canopy). The calculated diameter of the tower footprint for 
each site is listed in Table 3. To facilitate the comparison in this study, 
a cosine-law-based up-scaling method was applied to aggregate the 
30-m Landsat albedo to the tower footprint for individual sites (Shuai 
et al., 2011). The surface albedo corresponding to the tower footprint 



Fig. 6. Illustration of the aggregation from Landsat 30-m pixels (dotted gray grids) into the 
footprint projected on the ground by albedo meter (FOV = a) mounted on the tower h 
meters above canopy. 



Fig. 7. Comparison of actual (also called blue-sky) shortwave albedo between ground 
measurements (Y-axis) and satellite retrievals (X-axis) at 10:30 AM from Landsat re- 
trievals and at local solar noon from operational MCD43A3 (V005) products over six 
AmeriFlux forested sites in Pacific Northwest region of United States, both Landsat and 
MODIS meet the nominal 0.02 accuracy requirement in root mean square error (Sellers 
et al., 1995). The dashed lines represent an absolute accuracy of 0.03 compared to the 
ground data. 


ottower - footprint was obtained as the sum of available 30-m retrievals 
(ow(0) weighted by cos(0/), where 6 t is the view angle between the 
tower top and the center of pixel i (Fig. 6) for all the N pixels that fallen 
in the footprint of ground measurements (Eq. 12). 


_!£,(«(*) .a*©) 

Qt-tower -footprint v^N /n N 

Li=l COS ( 0 i) 


( 12 ) 


4.1.3. Comparison with ground measurement 

We derived the surface albedo from the 30-minute tower measured 
downwelling and upwelling radiation at 10:30 AM for Landsat and at 
local solar noon for MODIS over each site. Note that retrievals from 
Landsat as well as MODIS calculate intrinsic surface albedo under two ex- 
treme incident radiation situations (“black-sky albedo” corresponding to 
purely direct solar illumination and “white-sky albedo” corresponding 
to purely isotropic illumination), while the field measurements record 
the actual illumination corresponding to a mixture of both direct and dif- 
fuse radiation. To obtain comparable metrics with field measurements, 
we calculate the actual albedo (also called “blue-sky albedo”) via the in- 
terpolation between black-sky and white-sky albedos weighted by the 
ratio of direct or diffuse to the total downwelling radiation (Lucht et al., 
2000; Roman et al., 2011; Schaaf et al., 2002). Since the in-situ datasets 
lack information on direct/diffuse ratios, we simulated the direct/diffuse 
ratios for required solar zenith angles using 6S based on the simultaneous 
MODIS Terra atmosphere optical depths at the 550 nm band. Errors in- 
duced by the difference of defined wavelength interval for the shortwave 
band (ground 0.3-3.0 pm, Landsat 0.3-3.0 pm, and 0.3-5.0 pm for MODIS) 
are negligible because the solar irradiance beyond 2.5 pm accounts for less 
than 1.8% of the total between 0.3 and 14.3 pm (Hulstrom, Bird, & Riordan, 
1985). 

The scatter plot (Fig. 7) compares the Landsat blue-sky albedo aggre- 
gated to the tower field-of-view with the in-situ measured albedo at 
10:30 AM in the shortwave for the six AmeriFlux network sites. Re- 
trievals with snow and cloud contamination were removed from the 
analysis using the snow and cloud flags in the QA word of the satellite 
products. The Landsat retrievals are in very good agreement with the 
tower-based albedo, with a root mean square error (RMSE) less than 
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"pre-MODIS era" LUT approach L^d "Concurrent” approach 



Shortwave W$A from “pre-MODIS era” LUT approach 


Fig. 8. Cross-comparison of shortwave white-sky albedo map on day 2007-08-29 generated from “pre-MODIS era” LUT approach (left) with “concurrent” approach (middle), and the re- 
lated scatter plots over all available pixels (right). 


0.016 and a bias no more than 0.007. Discrepancy between the Landsat 
and ground albedos is confined to within ±0.03 albedo (dotted line in 
Fig. 7), which supports the absolute accuracy requirement (0.02-0.05) 
established by the climate modeling community (Sellers et al., 1995). 
Compared with the operational MODIS (V005) shortwave albedo re- 
trievals at local solar noon via the ground measurements as a bridge, 
the Landsat retrievals are slightly higher, except for the US-NR1 site. 
This makes sense if we consider the definition of black-sky albedo as de- 
scribed previously. Because values of black-sky albedo depend closely 
on the direction of solar illumination (i.e. solar zenith angle or timing 
of observation), and black-sky albedo is commonly observed to de- 
crease from sunrise to noon, then increase from noon to sunset, as val- 
idated for MODIS in Liu et al., 2009. 


4.2. Cross-comparison with the “ concurrent ” approach 

As an initial validation, we compared albedo maps generated by the 
“pre-MODIS era” LUT approach to those generated by the previously 
published “concurrent” approach, which has been validated previously 
(Roman et al., 2013; Shuai et al., 2011). Fig. 8 shows the shortwave 
broadband white-sky albedo maps derived from both approaches for 
the identical date 2007-08-29 (fill values excluded in the analysis). 
The “pre-MODIS era” LUT approach derieved albedo (left) is consistent 
with that from the “concurrent” approach (right) for the spatial varia- 
tion of albedo values, from low in the PNW forest region to the high in 
the eastern barren land. In general, albedo value extracted by the 
“pre-MODIS era” LUT approach is slightly higher than the “concurrent” 



Fig. 9. Illustration of the consistency between the “pre-MODIS era” LUT approach and the “concurrent” approach of Shuai et al. (2011). Shortwave white-sky albedo maps generated re- 
spectively from “pre-MODIS era” LUT approach (panel A) and “concurrent” approach (panel B), over an undisturbed forest region in Montana (fill value or disturbed regions in black). The 
absolute difference maps of white-sky albedo between “concurrent” and “pre-MODIS era” albedo (panel C, fill value or disturbed regions in white) for the overlapping years 2001-2011. 
Day of year is indicated for each albedo map (YYYY-DOY). 
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Fig. 10. Cross-comparison of shortwave White- Sky Albedo (WSA) map on day 2007-08-29 generated from “pre-MODIS era” LUT approach (left) with concident MCD43A3 (middle), and 
the related scatter plot over all available pixels (right). 


albedo in some scattered less homogeneous areas in areas of sparse veg- 
etation, which may be due to the difficulty in achieving high quality an- 
isotropy information in fairly heterogeneous localities. However, the 
pixel-value-based scatter plot (Fig. 8, right) exhibits the good agree- 
ment between both approaches with a high R 2 (0.988), a slightly posi- 
tive bias (0.005) and RMSE (0.009), including a small bias from 1-1 
line which is mainly induced by the lack of high quality anisotropy in- 
formation. The time series of shortwave white-sky albedo retrieved 
from both approaches at an undisturbed 30 km by 30 km forest region 
in Montana for the years 1985-2011 and 2001-2011 is displayed in 
Fig. 9A and B. The consistency over the overlapped summer days of 
2001-2011 is displayed in Fig. 9C. The albedo distribution map from 
both series varies simultaneously with the fluctuation of day-of-year 
(DOY), exhibiting lower values on day 2002-255 than other earlier 
DOY. The absolute albedo difference during 2001-2011 for this region 
is mostly less than 0.01. 


4.3. Cross-comparison with MODIS operational product 

To further examine the LUT performance, we compared an albedo 
map generated by the “pre-MODIS” BRDF LUT approach to that pro- 
duced by the operational MODIS albedo algorithm, which has been ex- 
tensively validated (Jin et al., 2003; Liu et al., 2009; Roman et al., 2009; 
Wang et al., 2014). To perform the comparison, the MODIS 500-m 
MCD43A (V005) product was reprojected to the Landsat UTM grid, 
and resampled to 30-m resolution. Pixels with fill values in either the 
MODIS or Landsat retrievals were excluded from the analysis. The oper- 
ational MODIS white-sky albedo for WRS-2 path/row 45/29 was pro- 
duced for day 2007-08-29 (DOY 241) using the instantaneous BRDF 
estimate inverted from the 16-day directional observations of days 
241-256. The resulting shortwave broadband white-sky albedo 
(Fig. 10) shows reasonable agreement between the “pre-MODIS era” 
LUT approach (left) and the MCD43A operational product (middle). 
The higher spatial resolution of the Landsat retrievals provides 
greater dynamic range in albedo values. For instance, scattered high 
values in upper right and central region may be compensated by the 
surrounding dominant lower values when aggregated to 500 m resolu- 
tion, resulting in the bright yellowish patch at 500-m resolution. In the 
scatter plot, the pixel-by-pixel comparison between the “pre-MODIS 
era” approach and MCD43A reveals the reasonable agreement with 
dominant points along the one-to-one line, but apparent variation due 
to the 30-m versus 500-m resolution. In addition, some variance may 
be introduced by the discrepancy of band passes between Landsat 
and MODIS. 


5. Discussion and conclusions 

Our previous study (Shuai et al., 2011) described a “concurrent” 
approach to retrieving 30-m albedo by combining Landsat directional 
reflectance with contemporary MODIS BRDF information. Here we 
have described an approach to extend the generation of albedo products 
back to the 1980s (“pre-MODIS era”) via construction of a priori anisot- 
ropy Look-Up Tables (based on modern MODIS data) that may be com- 
bined with historical Landsat reflectances. Preliminary LUTs for the 
Pacific Northwest (PNW) region of the United States were created, 
with BRDF information stratified by the combination of landscape attri- 
butes (land surface, terrain slope, season, and disturbance age and type) 
for the six non-thermal Landsat-5 and -7. These LUTs provide various 
BRDF trajectories for undisturbed landscape types and disturbed forest 
scenarios, including recovery trajectories for fire and non-fire disturbed 
forest regions. In the future, this approach for the creation of BRDF-LUT 
could be used to build up LUTs over other regions. With the assumption 
of generally invariant BRDF shapes for similar land cover conditions, the 
spectral white-sky and black-sky albedos are derived through albedo- 
to-nadir reflectance ratios as a bridge between the Landsat and MODIS 
scales, followed by a narrow-to-broadband conversion to produce the 
broadband albedos for visible, near infrared, and shortwave spectral 
regions. 

The accuracy of retrieved Landsat albedos were evaluated using 
available ground measurements at forested AmeriFlux stations in the 
PNW region and further evaluated through cross-comparison with the 
surface albedos generated from the “concurrent” approach as well as 
coincident MODIS operational products. The results show that the re- 
trieved Landsat 30-m shortwave albedo values meet the absolute accu- 
racy of 0.02-0.05 required by climate models at the validation stations 
(Sellers et al., 1995), and furthermore are consistent with those pro- 
duced from both “concurrent” approach and MCD43A algorithm. An ad- 
vantage of the Landsat albedo maps is that the finer resolution results 
include more detailed representation of the landscape, and a wider dy- 
namic range compared to the 500-m MCD43A product, especially in 
heterogeneous regions. This is particularly important given the fine- 
scale landscape patterns that result from both ecosystem disturbance 
and human modification of land cover. 

Future work will expand the validation to structural characteristics 
across the full range of the North American vegetation regimes by incor- 
porating recent Landsat-based disturbance products (Kennedy et al., 
2012; Masek et al., 2013). In addition, with the lengthening MODIS re- 
cord and the expansion of the detection of disturbance types in the re- 
cent disturbance products, we may be able to separate the non-fire 
disturbance BRDF-LUT into individual disturbance categories, and 
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reduce the need for BRDF interpolations. Equally important is the incor- 
poration of soil anisotropy once sufficient understanding of natural 
structured soil anisotropy features are obtained from not only laborato- 
ry measurements but also satellite-based estimates. Finally, we intend 
to develop an approach to retrieve snow-covered albedo at 30-m reso- 
lution to further understand the albedo evolution and energy budget 
over snow covered forest regions. 
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